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> 1 INTRODUCTION 



ABSTRACT 

In previous studies we identified two classes of starless cores, thermally subcritical 
and supercritical, distinguished by different dynamical behavior and internal structure. 
Here we study the evolution of the dynamically-unstable, thermally-supercritical cores 
by means of a numerical hydrodynamic simulation that includes radiative equilibrium 
and simple molecular chemistry. From an initial state as an unstable Bonnor-Ebert 
(BE) sphere, a contracting core evolves toward the configuration of a singular isother- 
mal sphere (SIS) by inside-out collapse. We follow the gas temperature and abundance 
of CO during the contraction. The temperature is predominantly determined by radia- 
tive equilibrium, but in the rapidly contracting center of the core, compressive heating 
raises the gas temperature by a few degrees over its value in static equilibrium. The 
time scale for the equilibration of CO depends on the gas density and is everywhere 
shorter than the dynamical timescale. The result is that the dynamics do not much 
affect the abundance of CO which is always close to that of a static sphere of the same 
density profile, and CO cannot be used as a chemical clock in starless cores. We use 
our non-LTE radiative transfer code MOLLIE to predict observable CO and N2H+ 
line spectra, including the non-LTE hyperfine ratios of N2H+, during the contraction. 
These are compared against observations of the starless core L1544. The comparison 
indicates that the dust in L1544 has an opacity consistent with ice-covered rather than 
bare grains, the cosmic ray ionization rate is about 1 x 10~^^ s~^, and the density 
structure of LI 544 is approximately that of a Bonnor-Ebert sphere with a maximum 
central density of 2 x 10^ cm~^, equivalent to an average density of 3 x 10^ cm~'^ within 
a radius of 500 AU. The observed CO line widths and intensities are reproduced if the 
CO desorption rate is about 30 times higher than the rate expected from cosmic-ray 
strikes alone, indicating that other desorption processes are also active. 



The starless cores are well described as small, dense, self- 
gravitating clouds supported largely but not entirely by 
thermal pressure. They are significant in the interstel- 
lar medium a s the future birthplaces of stars (reviews: 
iDi Francesco'et al. 2007: Bcrgin & Tafalla 20 j). In previ- 
ous papers (|Keto fc Fieldll2005l : lKeto fc Casellill2008l ) we de- 
veloped a simple model to describe the temperature, density, 
chemistry, and dynamical evolution of the starless cores. 
Those studies suggest that the starless cores can be con- 
ceptually divided into two classes, thermally subcritical and 
thermally supercritical depending on whether their central 
densities are very roughly less or greater than a few 10^ cm""^ 
respectively. This density is significant as the approximate 
value of two critical densities in the physics of starless cores. 
First, at this density gas cooling by coUisional coupling with 
dust is about as efficient as cooling by molecular line radia- 
tion. The gas temperature in the center of the supercritical 
cores is therefore significantly lower than in their envelopes 
and lower than found anywhere in the subcritical cores. Sec- 



ond, at this density, cores of a few Mq are at critical dynami- 
cal stability with respect to gravitational collapse. The ther- 
mally subcritical cores are stable against gravitational col- 
lapse and their dynamics may be d ominated by oscillations 
l|Lada et al.ll200a lAguti et al.ll2007l ) . In contrast, the super- 
critical cores are un stable with predominantly inw ard ve- 
locities (IWiUiams et al.l99g.; Casclli et al. 2002 : .Keto et all 
l2004l : ISohn et al.ll2007l : ISchnee et allbOOTf ) 

I n prev ious papers l|Keto et al.1 l2006l : iBroderick et al.l 

|2007| . |2008| ) we modeled the internal oscillations (sound 
waves) in thermally subcritical cores and computed some 
observable molecular spectral line profiles during the oscil- 
lations. In this paper, we model the dynamical and chemi- 
cal evolution of a thermally supercritical core and compute 
observable molecular spectral line emission during the early 
stages of gravitational contraction leading up to free-fall col- 
lapse. 

Our model is a one-dimen sional numerical hydrodynam- 
ics code (|Keto fc Fieldl I2OO5I ) that includes the radiative 
equilibrium of dust, gas cooling by both molecular lines and 
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coUisional cou plinR with the dust, and a simphfied molecu- 
lar chemistry (|Keto fc Caselli|[2008l ). To compare the model 
with observations we use our numerical code for non-LTE 
radiative transfer MOLLIE t o predict th e observable molec- 
ular spectral line emission (|Ketol Il990l : iKeto et alJliooi ). 
Previously we modeled overlapping hyperfine emission by 
assuming that the hyperfine components are in statistical 
equilibrium. Here we use an updated algorithm that does not 
make this assumption and better reproduces the non-LTE 
hyperfine line ratios (excitation anomalies) th at are seen in 
N2H" ' " spectra from man y of the starless cores l|Caselli et alj 
1 19951 : iDaniel et al.ll2007l ). We compare these predictions to 
previous observations of one specific core, L1544 whose prop- 
erties place it in the thermally supercritical class. 

There has been some considerable research into the 
topic of how gas clouds contract to form pro tostars inclu ding 
some classic papers in the fiel d (|Henvev et a l. 1955: Ha vashil 
Il96ll : iLarsoiJl 19691 : [siuj|l977l '). More recent studies also ad- 
dress the chemical evolution of the starless cores using a 
variety of models for the dynamics of gravitationally con- 
tracting cores^ 

.Aikawa et al.l (|200ll . 120031 ) followed the chemical evolu- 
tion of a gravitationally collapsing core assuming the Larson- 
Penston (LP) approxi mation for the densit ies and velocities. 
In a subsequent paper. [Aikawa et al. I (j2005) imp roved on this 
mode l using a numerical hydrodynamic code (Ogino ct al. 
Il99g| ) to follo w the collapse of an isothe rmal Bonno r- Ebert 
(BE) sphere. lEvans et all l|200ll ') and iGalli et all (|2002l ') 
showed that in the case of static BE spheres there was little 
difference between the structure of an isothermal BE sphere 
and one in radiative equilibrium. Our numerical hydrody- 
namic simulation that allows for radiative equilibrium dur- 
ing coll apse also confirms th at the isothermal approximation 
used in lAikawa et al.l (|2005l ) is quite adequate for dynamic 
as well as static BE spheres (our i]2.2|) . However, if we go on 
to predict molecular line strengths from the model cores, we 
need to consider temperature variations. Our discusses 
the spectral line model i ng. 

iRawlings fc Yate^ (|200ll ) used the similarity solution 
of th e collapse of a singular isothermal sphe re (SIS) (IShul 
I1977I) in their study of chemical evolution. iTsamis et al.l 
l|2008t ) combined this model for the dynamics with a model 
for the internal temperature structu re of a BE sphere from 
lEvans et al.l (120011 ) . iLee et ahl (|2004l ) approximated the con- 
traction of a BE sphere as a progression of static spheres 
of increasing central density. Our hydrodynamic simulation 
suggests that this series of static equilibria is a good approx- 
imation during the subsonic phase of contraction, and that 
the SIS model is a good approximation at the endpoint of 
this phase (our 

IVan Weeren et al.l (|2009l ) and lBrinch et al.l (|2008l ) mod- 
eled the chemical evolution during the formation of a 2D 
rotating accretion disk and protostar. Our spherical model 
applies to earlier evolutionary times while the contracting 
core still maintains its B E structure. 

iLi et al.l l|2002l ) and IShematovich et all (|2003l ) followed 
the chemical evolution during the contraction of cores that 
are supported by a spherical approximation of magnetic 
pressure. The contraction is controlled by leakage of the 
magnetic field out of the core by ambipolar difi'usion. In 
the spherical approximation, this is modeled by a diffusion 
equation in Lagrangian coordinates (equation 3 of iLi et al.l 



|2002|) . One difficulty with assessing the applicability of 
spherical models of magnetic collapse is that the approxi- 
mation of spherical symmetry results in a density structure 
that, similar to the SIS and BE spheres, can be approxi- 
mated by a polytrope with an index that depends on the 
effective equation of state. Observationally, it is difficult to 
distinguish polytropes of different index because their den- 
sity structures are most different from one another at small 
and large radii where the observations have difficulties with 
angular resolution and sensitivity (low density gas) respec- 
tively. 

In contrast, there is strong observational motivation 
for the model of thermal support. Masses, densities, and 
temperatures estimated in surveys of starless cores sug- 
gest that thermal energy provides a bout 75% of the in- 
ternal energy required for equilibrium (Pic kman fc Clemens! 
|l983; Myers fc Benson 1983; Lada ct al. 200j). The remain- 
ing 25% could be in magnetic energy or wave energy or a 
combination of both. While the 25% non-thermal energy is 
important for the dynamical stability and evolution of the 
core, the structure of the core is determined largely by the 
dominant thermal energy. 

Our study is different from previous studies in that we 
consider the evolution of a BE sphere rather than an accre- 
tion disk plus protostar, or an SIS sphere, or a magnetically 
supported sphere. Our simulations is the first that follows 
the dynamics of the contraction simultaneously and self- 
consistently with radiative equilibrium and a simple model 
for the chemistry that allows for changes in the molecular 
line cooling as the molecular abundances in the gas phase 
change with time. 

We first discuss the dynamics of the gravitational con- 
traction of a BE sphere. Our code is best suited for subsonic 
velocities and we confine our discussion to the early phase 
of evolution. We then discuss how the gas phase abundance 
of CO evolves during contraction. We predict how the 
molecular line strengths should vary across the model from 
center to edge at different times during the contraction, 
and we compare this prediction with ratios of line strengths 
previously observed in L1544. We intend the comparison 
in a general sense only and not in detail. For example, the 
model is spher ically symmetric while both the ob served 
morphologv (IWard-Thompson. Motte. fc Andrd Il999l: 
iBacmann ct al."2000') and observed spectra ("William s et al.l 
1 1999 ; Casolh ct al. 2002; Willi ams et a l. 2006) indicate that 
both the density and velocity structure of L1544 are more 
complex than can be reproduced by a purely spherical 
model. 



2 THE GRAVITATIONAL CONTRACTION OF 
AN UNSTABLE BONNOR-EBERT SPHERE 

We model a core of 10 Mq that begins its evolution with an 
initial density of 2 x 10* cm~^ in a state of radiative and dy- 
namical equilibrium. The mass and density place the core in 
the thermally supercritical class (figure 14 o f iKeto fc Caseiiil 
12008) so that the dynamical equi l ibrium is unstable. Based 
on the results of iKeto fc CasellH (|2008l ) we use our higher 
dust opacities t hat are equal to four t imes the standard 
dust opacities of lOssenkopf fc Hennind (|l994l ). and we use 
a rate of cosmic ray ionization of 1.3 x lO"^'^ s~^ referred 
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to as the "low" rate in lKeto fc Casellil (120081 ). In that inves- 
tigation, this combination was found to result in gas tem- 
peratur es that match those suggested by the observations 
of .Crapsi et alj (|2007l ) at their suggested central density of 
2 X 10^ cm"^. 

Figures [1] and [2] show the density and velocity profiles 
during contraction as the central density increases from its 
initial value to 2 x 10* cm~^ over a time period of 1 Myr. The 
shapes of the density and velocity profiles that are generated 
by the numerical simulation can be understood in terms of 
relatively simple theory. We begin with the density profile. 



2.1 The density profile 

The density profiles in figure [l] show an outer region 
where the density scales as and an inner region 

where the density is nearly constant. As the evolving core 
contracts, the density profile maintains this shape even as 
the central density increases (^Larson '1969*; 'Hunt eij|l977l: 
Foster fc Chevalier 1993; Ogin o et al.. .1999 : Kan dori et al.l 



20051 ). Thus we should expect to observe this characteristic 
shape in almost all starles s cores, and this is g e nerall y 



i g e nerall y 

confirmed b y observ atio ns. IWard-Thompson et all Il993), 
lAndre et al] (llQQd). IWard- T hompson . Motte. fc Andrei 
and iBacmann et al. suggested that the 



density profiles derived from observations of dust emission 
could be well matched by two power laws, one shallow and 
one steep. This approximation is also suggested by the 
numerical st udies of the c ontraction of a BE sphere cited 
above. Tafalla et al.1 (|2002| . §4.2) suggested that the density 
profile can be better fit by. 



p{r) = 



(1) 



1 + (r/rf)- 

where r/ is the radius of the inner region. This equation 
along with r/, which we define below, provides a convenient 
approximation for future use. 



2.1.1 The flat inner region 

In the center of a BE sphere, where the sound-crossing 
time is less than the free-fall time, density perturbations are 
rapidly smoothed by pressure waves. Also the self-gravity in 
the center is relatively weak so that the gas is essentially 
confined by the pressure of the overlying gas. Both these 
effects act to maintain constant density around the center. 

Previously, the size of the region with a flat density 
profile, r/, was determined empirically for individual cores 
from their observations. We can also determine this radius 
theoretically. At r/ the ratio of the sound-crossing and free- 
fall times ought to be one. Therefore, r/ is the product of the 
sound speed, a, and the free-fall time at the central density. 



(32p,G/37r)i/2- 



(2) 



This provides a better approximation in equation 1 than 
the standard s cale length of a BE sphere, [kT/ (AnGpc)]'''''^ 
llBonnoJ 1195^ 1 and is also useful in suggesting the physics 
behind the behavior of a contracting BE sphere. The ra- 
dius, rf is equal to the Jeans length for sound speed, a, and 
density, pc. 



Figure \T\ compares the density profile computed by 
equations [T] and (2] assuming a gas temperature of 10 K and 
a central density of 2 x 10^ cm~'^ (green curve), against the 
density profile computed by the numerical hydrodynamic 
simulation. The comparison suggests that equations [1] and [2] 
provide a good, simple approximation to the density profile 
of a BE sphere. (A more accu r ate ap proximation is described 
in lNataraian &: Lvnden-Belll (|l997l ).) Because the numerical 
hydrodynamic simulation also accounts for a variable tem- 
perature set by radiative equilibrium, the comparison in fig- 
ure 1 also shows that the departures from isothermal owing 
to radiative equilibrium have little effect on the overall den- 
sity structure of the core. 

Figure [1] also shows the characteristic radius, r/, for 
the different evolutionary times during the contraction. The 
figure shows that equation [2] provides a good approxima- 
tion to the turn-over radius, r/, at all times in the evo- 
lution. According to equation 2, this width shrinks as the 
central density of a contracting BE sphere increases. In the 
subsonic phase of gravitational contraction, r/ moves at 
about half the maximum gas velocity where this maximum 
is taken from the velocity profiles for each evolutionary time 
as shown in flgure[2] 

In comparing our theoretical results to several observa- 
tions, we find that a central density of about 2 x 10^ cm"'' 
provides a better match to the observations (f|4]) than cen- 
tral densities that are a factor of 10 higher or lower. This 
estimate is in agreement with th e density of 2 x 10^ cm"'' 
suggested by the observations of ICrapsi et al.l l|2007l ) if we 
average the theoretical density profile over the size of their 
observing beam. Although the theoretical density profile ap- 
pears flat in the center of a log-log plot of density versus 
radius, it is sharply peaked within the size of the observing 
beam. The average density within a radius of 500 AU (ob- 
serving beam size of 1000 AU = 7" at 140 pc) is 2.8 x 10*^ 
corresponding to a peak density of 2 x 10*^ cm"''. In fact, at 
a spatial resolution of 1000 AU, observations have difficulty 
measuring the central density. For example, the average den- 
sity over 1000 AU corresponding to a peak density of 2 x 10* 
cm"'' is only 3.4 x 10® cm"'' , little different from the average 
density of 2.8 x 10® cm~^ corresponding to a peak density 
of2xl0^cm"^ 

In our previous paper (|Keto Sz Casellil I2OO8I ). we used 
the lower central density of 2 x 10® cm"^ and stated that a 
higher central density would imply a higher density through- 
out the core. This is not correct. Owing to the inside-out 
character of the collapse, the central density increases much 
more rapidly than the densities outside of r/. As explained 
in the next section, in an evolving core, the densities over 
most of the core are actually about the same for both these 
higher and lower central densities (figure [T|. 



2.1.2 The self-similar outer region 

In their numerical hydrody namic simulations of gravi t ation - 
ally collapsing gas clouds, iBodenheimer fc SweigartI (|l968h 
found that regardless of the initial configuration of the cloud 
or the initial conditions, all their examples evol ved to den - 
sity profiles scaling as r~^ in their outer regions. IShul (|l977l ) 
suggests that this scaling is a property shared by all self- 
gravitating isothermal systems under certain general condi- 
tions. The system should evolve subsonically, meaning that 
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it should be close to hydrostatic equilibrium initially, and 
the outer region should not be affected by an outer bound- 
ing surface. These conditions are applicable to contracting 
BE spheres. The initial state is one of hydrostatic balance, 
albeit unstable. Although BE spheres have an outer bound- 
ary which is the radius at which the Lane-Emden equation 
is truncated, the external pressure at the boundary is the 
same as would be provided by the Lane-Emden equation if 
the solution were continued. Thus in the initial configura- 
tion, the density profile is the same as if the core had no 
outer boundary. 

Figure [2] shows that the BE sphere evolves by subsonic 
contraction to resemble the SIS, a result also foun d in ear- 
lier studies (|Larsonlll969l : [Foster fc Chevalieilll993l '). As the 
central density increases with the contraction, the character- 
istic radius, r/, moves inward, the constant density region 
shrinks in size, and the outer region with its density profile 
scaling as p ~ r~'^ accounts for more and more of the core. 
Because the density of the outer region changes very slowly, 
the density contrast between the inner and outer regions also 
increases rapidly. Thus, as the contraction proceeds, the BE 
sphere evolves to resemble an SIS with its density profile 
everywhere and its infinite density contrast. 

This behavior is expected if we recall that the static 
equilibrium solutions of the Lane-Emden equation form a 
continuous series with increasing density concentra tion an d 
the SIS as the limiting case (,Chandrasekhar..l957. ; ShuligrTh . 
During the subsonic, quasi-static phase of contraction, an 
unstable BE sphere evolves approximately along this series 
of static equilibria toward the limiting case of an SIS. 



2.2 The velocity profile and the inside-out 
collapse of a BE sphere 

The velocity profile during the subsonic phase of the con- 
traction of a BE sphere (figure 2) shows a characteristic 
A-shape with the inward velocity as a function of radius in- 
creasing from near zero in the core center to a maximum 
just outside the characteristic radius, rf, before decreas- 
ing again toward the edge of the core (iFoster fc Chevalierl 
ll993l : IOgino et al.lll999l : iKandori et al.ll2005l ). The origin of 
this profile can be understood by considering the forces in- 
side and outside of r/. In the very center, the velocity is 
zero by symmetry. In the region, rf, where the density is 
constant, the instantaneous acceleration due to the gravi- 
tational force increases linearly with radius. The velocity, 
which is the time-integrated acceleration, may also be ex- 
pected to increase with radius as well, and this is verified 
by the numerical evolution. Outside of r/, the density falls 
off as . Here the instantaneous gravitational accelera- 
tion decreases with radius along with its infiuence on the 
velocity. During the contraction, the redistribution of mass 
in the center does not change the gravitational acceleration 
in the outer region because of spherical symmetry. Thus, in 
the outer region, the only change to the hydrostatic balance 
is due to the pressure force which propagates outward from 
the center only at the sound speed. These several effects that 
would change the hydrostatic balance in the outer region all 
decrease as we move outward, away from rf. Thus the outer 
region remains in approximate hydrostatic balance with low 
inward velocities for about a sound-crossing time, on the or- 



der of 1 Myr, and the velocities are lower at larger radii, 
further from r/. 

The same characteristic A-shaped profile also develops if 
we consider the pressure-free collapse of a BE sphere. In this 
case the gas pressure is ignored and the velocity field is deter- 
mined s olely by the variation of the gravitational force wit h 
radius. (jWhitworth fc Ward-Thompsonll200ll : lMversll2005l ). 
Of course, the gas pressure is required to obtain the initial 
BE density distribution in the first place. 



3 THE TEMPERATURE AND CHEMISTRY 
OF THE EVOLVING CORE 

In the previous section we saw that the dynamics of a con- 
tracting BE sphere, as followed by a numerical simulation 
that included radiative and chemical equilibrium, were well 
matched by a simpler isothermal approximation. In this sec- 
tion we will see that the temperature and chemistry of 
an evolving core have a significant effect on the observ- 
able molecular line spectra. The gas temperature directly 
affects the line brightness through the source function. The 
chemistry affects the line brightness through the gas phase 
abundance of the observable molecules. The most signifi- 
cant processes affecting the abundances are the depletion 
of molecules from the gas phase as they freeze onto dust 
grains in the dense center of a core and photodissociation 
of molecules near the core boundary. We calculate the CO 
abund ance with a simple mo del that includes these two pro- 
cesses (|Keto fc Casellill2008l ). 

Figure |3] shows the gas temperatures, density, velocity, 
and CO abundance at 3 different times during the collapse 
when the central density is 2 x 10^, 2 x 10^, and 2 x 10* cm~^. 
As discussed further in i|4] , a model with a central density 
of 2 x 10^ cm~^ provides a better match to the observations 
than central densities that are a factor of 10 lower or higher. 
Also shown are the density, temperature, and abundance 
of static cores that have the same central densities as the 
dynamically evolving core. The properties of the contracting 
cores are similar to those of the static cores with the same 
central density, with the exception of the velocity. 



3.1 The temperature in the center of a 
contracting core 

The starless cores are heated from the outside by both cos- 
mic rays and by the diffuse interstellar radiation field sup- 
plied by the population of stars in the Galaxy. The cores 
are cooled by the long wavelength radiation from cold dust 
and by molecular lines. Because the incoming shorter wave- 
length starlight is strongly absorbed by the cores whereas 
the outgoing long wavelength radiation from cold dust is 
not, the denser cores are coldest in their centers. This tem- 
perature structure has previously been calculated in q uite 
a number of papers (references in iKeto fc Fieldl |2005| ). In 
figure [3] we see that the central temperature of the static 
core with a central density of 2 x 10** cm~^ is quite low, 
below 5 K. In the evolving core, compressive heating keeps 
the temperature above 5 K, more consistent with the tem- 
peratures inferred from molecular line a nd dust observations 
(jCrapsi et al.ll2007l : [Pagani et al.ll2007l ). 
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3.2 CO depletion and desorption 

In our previous paper we assumed tliat tlie rate for tiie 
desorption of CO off dust, tlie inverse process to deple- 
tion or freeze-out, was due to whole-grain iieating b y cosmic 
rays (equation ll lKeto fc Cas clli 2008; H asegawa fc HerbstI 
1 19931 ). When we compare the C^*0 and C^^O spectra pre- 
dicted by our radiative transfer code MOLLIE with observa- 
tions of L1544 we find that the predicted lines are not bright 
enough to match the observations. Since these CO lines are 
optically thin, the CO column density is too low. 

The desorption rate cannot be increased by simply in- 
creasing the flux of cosmic rays. In our earlier study we 
found that a higher flux of cosmic rays would cause the gas 
temperature at mid-r adii to be higher tha n suggested by 
observations (figure 12 lKeto fc Casellill2008l ). 

There are other processes in addition to direct cosmic- 
ray strikes that caus e desorption and increase t he gas phase 
aburi d ance of CO JPrasad fc Tarafdai _ 19831: Leger et ahl 



JCaselli et al.l I2OO2I: iKeto et al.l [iooi: iBergin et al.l I2OO6I : 



19851: Id'Hendecourt et a.lj Il985l: iDulev fc Williamj Il993l : 
Willacy fc Millar' Il998l: iTakaha shi fc Williams! I2OO0I: 



Dalg arno 2006; O berg et ah 2009). S hen et all l|2004l ) found 



that the energy transferred to dust by the UV radiation 
field produced by cosmic-ray strikes on molecular hydrogen 
is almost one order of magnitude larger than the energy 
transferred directly to d ust by the cosmic-ray particles 
themselves. iRoberts et al.l |2003) suggested that in addition 
to direct and indirect heating by cosmic rays, another 
significant source of heat is the exothermic formation of H2 
on the grain surfaces. [Roberts et all (|2007h suggest that the 
rates of these heating processes are not known. However, 
because all the processes depend on the first power of the 
density, the same as direct cosmic ray heating, we can 
account for additional desorption processes in our model by 
simply increasi ng the desorption rate above that given by 
equation 11 of iKeto fc Casellil (120081 ) while stiU assummg 
the same first-power dependence on the gas density. 

Figure 3] shows the abundance of CO obtained as an 
equilibrium between the 4 processes of depletion, desorp- 
tion, photodissociation and chemical formation, and calcu- 
lated for 4 different desorption rates. The lowest rate is 
equivalent to desorption ca used by direct cosmic-ray strikes 
(|Hasegawa fc HerbstI \l99^ on dust at the "low" cosmic- 
ray rate (1.3 x 10~^ ^ s~^) as defined in equation 11 of 
iKeto fc Casellil (|2008h . The higher desorption rates are fac- 
tors of 3, 9, and 30 higher than the lowest rate. In order to 
match the observed CO spectra (figure [9]) we find that we 
need to increase our desorption rate by a factor of 30. At 
this rate, desorption and depletion have equal time scales at 
a density of about 10* cm~^. 

As shown in figure |4l with the highest desorption rate, 
the maximum CO abundance is about a factor of 2.5 higher 
than with the lowest rate. This increase in abundance makes 
the optically thin C^^O and C^*0 lines almost a factor of 2 
brighter and a better match to the observations. 

The gas phase abundance could also be increased by 
decreasing the reverse reaction of freeze-out. One way 
this could occur is if the dust grains in the center of 
L1544 were fewer in number but larger in size, for ex- 
ample by coagulation. This wou ld decrease the total sur - 
face area available for freeze-out (Ivan Dishoeck et~al]| 19931 ). 
There is some independent evidence for coagulated grains 



iFlower et al.H2006l : IVan Weeren et al.ll2009l ). However, to re- 
duce the rate of freeze-out by a factor of 30, would require 
a mean grain diameter of about 10 fim, too large accord- 
ing to available models of du st coagulation in dense clouds 
(|Ossenkopf fc Henning||l99^. Decreased depleti on at a rate 
slower than equation 8 of Keto fc Casellil (|2008l ) could con- 
tribute to the higher CO abundance, but could not bring 
the CO abundance up high enough to match the spectral 
line observations. Finally, decreased photodissociation might 
also increase the CO abundance. This could come about if 
the UV fiux were diminished, for example, if the core were 
embedded in a larger scale molecular cloud. We ran a sep- 
arate radiative transfer calculation and found that the CO 
abundance cannot be increased deep enough into the core 
to make a difference to the line brightness. 



3.3 Depletion time scale 

T he time scale for CO to freeze onto dust. Ton (equation 
8 iKeto fc Casellil I2OO8I : iRawlings et~aLlll992f) and the time 
scale for the inverse process of desorp tion, Tpff (30x equa- 
tion 11 Keto fc Casolh 2008; Hasega wa fc Herbst 1993) may 
be combined to deter mine the time scal e for the change of 
CO in the gas phase (|CaseUi et al.ll2002l ). 



TCO = 



TonTpff 
Ton + Toff 



(3) 



This time scale for the equilibration of the CO abun- 
dance varies across the core but is everywhere aster than 
the dynamical time. In the core center, the free-fall time is 
a close approximation whereas in the more slowly evolving 
outer region, the sound crossing time is more appropriate. A 
starless core with a mass of 10 Mq and a central density of 
2 X 10^ cm~^ has a free-fall time, t// = 0.03 Myr using the 
central density in the standard equation whereas the sound 
crossing time is about 2 Myr. In comparison, at the center of 
the core, the CO equilibration time scale tco ~ 0.01 Myr, 
a factor of 3 shorter than the free-fall time. At the edge of 
the core where the density is 20 cm~^, tco ~ 0.5 Myr, still 
a factor of 4 shorter than the sound crossing time, through- 
out the contraction. Figure [3] shows that the abundance of 
CO during the contraction of a BE sphere is everywhere 
quite close to its steady state value. For this reason, it is not 
possible to use CO as a chemical clock in starless cores. 



3.4 CO abundance and time reversibility 

The previous discussion demonstrates the importance of des- 
orption in determining the CO abundance. The desorption 
rate is required to determine the equilibrium abundance to- 
ward which the system evolves asymptotically. In particu- 
lar, the CO abundance should not be calculated simply from 
the depletion rate as an exponentially decreasing function of 
time as would be the case if desorption were ignored. Even in 
gas dense enough that the depletion rate is much faster, gen- 
eral considerations always require both the forward and re- 
verse reactions. Real physical systems must satisfy detailed 
balance and time reversibility (|OnsageI^ll931^ . A system de- 
scribed by a single rate equation, for example depletion only 
without desorption, does not satisfy these conditions, cannot 
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describe a physical system, and this description wiU produce 
misleading results. 



4 COMPARISON WITH OBSERVATIONS 

All other things being equal, a collisionally excited, optically 
thin line should increase in brightness with the path length 
and the gas density. Thus molecular lines should be bright- 
est through the center of the core. Observations of starless 
cores do not always show this behavio r because of varyin g 
molecular abundances within the core dXafalla et al.ll2002l l. 
For example, figure 2 of lCaselli et al.l (|l999l ) reproduced here 
as figure [S] shows the brightness of the (1-0) transitions of 
C^''0, C^*0, and N2H+ at a set of positions across the L1544 
core from edge to center to edge. The N2H''' molecule does 
not suffer much, if any, depletion from the gas phase at high 
densities, and the observed N2H^ brightness increases to- 
ward the core center whereas the observed CO brightness 
does not. This difference is due to the near complete deple- 
tion of CO at higher density, ft is possible th at NgH"*" may 
be sl ightly depleted in the center of LI 544 jCaselli et al.l 
|2002| ) although we do not include this effect in our model. 
We will investigate this possibility in a future modeling of 
higher density tracers and transitions. Our modeling here 
shows that the decrease in temperature in the core center is 
not sufficient to significantly decrease the line brightness. 

Figure |6] shows the simulated integrated spectral line 
brightness as a trace of positions across the model evolving 
cores in the same format as the observational data in figure 
[5l The integrated intensities include the emission from all 
three hyperfine lines of C^'^O(l-O) and all seven hyperfine 
lines of N2H^(l-0). At central densities greater than 10® 
cm~'^ the modeled brightness begins to look like the obser- 
vations. There are some difi'erences. 

The observational data might show a slight dip in the 
integrated intensity of the CO lines toward the center which 
is not present in the model. However, there is a decrease in 
the peak intensity of the modeled CO lines (figure [T}. This 
dip does not show up in the integrated intensity because the 
line widths in the model increase toward the center owing to 
the increase in the inward velocities in the inside-out collapse 
(figures |8]). The increase in the line width of CO is less than 
that of N2H^ because, owing to depletion, there is essentially 
no gas phase CO in the center of the core where the velocities 
are highest. From a density of 2 x 10® to 2 x 10** cm~^ the 
C"0 and C^*0 line widths increase by 0.07 and 0.06 kms~^ 
while the N2H+ line width increases by 0.15 kms~^. 

Because the density, temperature, and abundance struc- 
tures of contracting cores are so similar to those of static 
cores, the primary observable evidence of contraction is not 
the line brightness, but rather the shapes of those spectral 
lines such as N2H'''(l-0) that are optically thick enough to 
show self-absorption. Figure[5]shows the spectral line profiles 
of the (1-0) transitions of C^^O, C^**0, and N2H+ expected 
from our model for the 3 evolutionary times corresponding 
to figure [31 As the inward velocities increase during the con- 
traction, the N2H''' spectral lines become wider (figure [S]), 
and eventually split into two peaks (figure [9J because of the 
inward velocities in the core (figure [3}. In contrast, because 
there is little CO in the center where the velocities are high- 
est, the width of the CO lines changes very little (figure[8]) as 



the core evolves, and the CO lines do not split (figure[9} . The 
brightness of both the N2H"'" and CO lines are nearly con- 
stant as the central density increases from 2 x 10® to 2 x 10* 
cm^"* because there is very little change in the density of 
most of the core except for the center where CO is depleted. 
The slight decrease in the peak brightness of C^^O (figure 
[7] and [9]) is consistent with a constant integrated intensity 
and the slight increase in line width (figure [8} . As the core 
evolves from a central density of 2 x lO'^ to 2 x 10* cm~'^ 
there is very little change in any of the observed spectra, ei- 
ther CO or N2H^, because the density increase is happening 
in a very small region in the center that does not contain 
much mass relative to the entire core and is also becom- 
ing small compared to the observing beam. Thus the N2H"'" 
spectra (figure [2| for the last two evolutionary times look 
almost the same. Interferometric observations of higher den- 
sity transitions of N2IT*' and N2D^ will be needed to unveil 
the dynamical stage, structure and kinematics of the rapidly 
contracting center of the core. This will be investigated in a 
future paper. 



5 CONCLUSIONS 

We follow the contraction of a thermally supercritical core 
through the evolutionary phase of subsonic contraction with 
a numerical code that includes radiative and chemical equi- 
librium. 

We identify a characteristic radius, the product of the 
sound speed and the free-fall time, as the point where the 
density profile of a BE sphere transitions from an outer re- 
gion where the density scales as r~'^ to an inner region of 
constant or fiat density. In the inner region, the dynamical 
time is the gravitational free-fall time at the central density. 
In the outer region, the dynamical time is the much slower 
sound- crossing time. The difference results in inside-out col- 
lapse. 

Once this characteristic radius becomes smaller than 
the angular resolution, observations have difficulty in de- 
termining the exact central density even if the observed 
molecule is not significantly depleted. 

We follow the gas temperature and abundance of CO 
during contraction. In the rapidly contracting center of the 
core, compressive heating raises the gas temperature by a 
few degrees over its value in static equilibrium and the re- 
sults are consistent with observations. The temperature vari- 
ations from the radiative equilibrium do not significantly 
affect the dynamics of the large scale structure of the ther- 
mally supercritical core. The time scale for the equilibration 
of CO is everywhere shorter than the dynamical timescale. 
Thus, the CO abundance is always close to that of a static 
sphere of the same density profile. Therefore, the CO abun- 
dance cannot be used a chemical clock to estimate the age 
of a starless core. 

The comparison with line profiles observed toward 
L1544 suggests that the dust has an opacity consistent 
with ice covered grains, the cosmic ray ionization rate is 
close to lxlO~^^ s~^, and the maximum density is about 
2x10^ cm~^, equivalent to an average density of 3x10® 
cm""* within a radius of 500 AU corresponding to an ob- 
serving beam size of 1000 AU. The line width and intensity 
of C^*0 and C^'^O lines can be simultaneously reproduced 
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by our model only if the CO desorption rate is about 30 
times larger than the rate expected from cosmic-ray impul- 
sive heating alone. 



Supercritical Starless Cores 7 



© 0000 RAS, MNRAS 000, 000-000 



8 Keto and Caselli 



REFERENCES 

Aikawa, Y., Ohashi, N., Inutsuka, S.-L, Herbst, E., 

Takakuwa, S., 2001, ApJ, 552, 639 
Aikawa, Y., Ohashi, N., & Herbst, 2003, ApJ, 593, 906 
Aikawa, Y., Herbst, E., Roberts, H. & CaseUi, P., 2005, 

ApJ, 620, 330 

Aguti, E.D., Lada, C.J., Bergin, E.A., Alves, J.F., Birkin- 

shaw, M., 2007, ApJ , 665, 457 
Andre, P., Ward- Thompson, D., Motte, P., 1996, AA, 314, 

625 

Bacmann, A., Andre, A. P., Puget, J.-L., Abergel, A., Bon- 
temps, S., & Ward-Thompson, D. 2000, A&A, 361, 555 

Bergin, E. A., Maret, S., van der Tak, F. F. S., Alves, J., 
Carmody, S. M., & Lada, C. J. 2006, ApJ, 645, 369 

Bergin, E.A. & Tafalla, M., 2007, ARAA, 45, 339 

Bodenheimer, P. & Sweigart, A., 1968, ApJ, 152, 515 

Bonnor, W., 1956, MNRAS, 116, 351 

Brinch, C, Van Weeren, R.J., & Hogerheijde, M.R., 2009, 
AA, 489, 617 

Broderick, A. E., Keto, E., Lada, C. J., & Narayan, R., 

2007, ApJ, 671,1832 

Broderick, A. E., Narayan, R, Keto, E., & Lada, C. J., 

2008, arXiv.0804.1790 

Caselh, P., Myers, P., Thaddeus, P., 1995, ApJL, 455, L77 
Caselh, P., Walmsley, C, Tafalla, M., Dore, L., & Myers, 

P. 1999, Ap Jl, 523, 165 
Caselli, P., Walmsley, C, Zucconi, A., Tafalla, M., Dore, 

L., & Myers, P. 2002, Ap J, 565, 344 
Chandrasekhar, S., 1957, Stellar Structure (New York: 

Dover Publications) 
Dickman, R.L., & Clemens, P., R.L., 1983, ApJ, 271, 143 
Crapsi, A., Caselh, P., Walmsley, M. C, & Tafalla, M. 2007, 

AA, 470, 221 
Dalgarno, A., 2006, PNAS, 103, 12269 
Daniel, F., Cernicharo, J., Roeff, E., Gerin, M., Dubernet, 

M.L., 2007, ApJ, 667, 980 
D'Hendecourt, L. B., AUamandola, L. J., & Greenberg, 

J. M. 1985, AA, 152, 130 
Di Francesco, J., Evans, N.J., Caselli, P., Myers, P.C., 

Shirley, Y., Aikawa, Y., Tafalla, M., 2007, in Protostars 

and Plabets V, eds. B. Reipurth & D. Jewitt, Univ. of 

Arizona Press, Tucson, page 17 
Duley, W. W., & Wilhams, D. A. 1993, MNRAS, 260, 37 
Evans, N., Rawlings J., Shirley, Y., & Mundy, L. 2001, Ap 

J, 557, 193 

Foster, P. & Chevalier, R., 1993, ApJ, 416, 303 

Flower, D. R., Pineau Des Forets, G., & Walmsley, C. M. 

2006, AA, 456, 215 
Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, 

ApJ, 262, 590 

Galh, D., Walmsley, M., & Gongalves, J. 2002, AA, 394, 
275 

Hasegawa, T., & Herbst, E., 1993, MNRAS, 261, 83 
Hennebelle, P., 2003, AA , 411, 9 
Hayashi, C, 1961, PASJ, 13, 450 

Henyey, L.G., Lelevier, R., & Levee, R.D., PASP, 67,154 
Hunter, C, 1977, ApJ, 218, 834 

Kandori, R., Nakajima, Y., Tatematsu, K., Aikawa, Y., 
Naoi, T., Sugitani, K., Nakaya, H., Nagayama, T., Na- 
gata, T., Kurita, M., Nagashima, C, Sato, S., 2005, AJ, 
130, 2166 



Keto, E., 1990, ApJ, 355, 190 

Keto, E., Rybicki, G., Bergin, E., Plume, R., 2004, Ap J, 
613, 355 

Keto, E. & Caselh, P., 2008, arXiv0804.0822 
Keto, E. & Field, G., 2005, ApJ, 635, 1151 
Keto, E., Broderick, A., Lada, C.J., & Narayan, R., 2006, 
ApJ, 652,1366 

Lada, C. J., Bergin, E. A., Alves, J. F., & Huard, T. L. 

2003, ApJ, 586, 286 
Lada, C.J., Muench, A. A., Rathborne, J., Alves, J.F., & 

Lombardi, M., 2008, ApJ, 672, 410 
Larson, R.,B., 1969, MNRAS, 145, 271 
Lee, J.-E., Bergin, E.A., Evans, N.J., 2004, ApJ, 617, 360 
Leger, A., Jura, M., & Omont, A. 1985, AA, 144, 147 
Li, Z.-Y., Shematovich, V.I., Wiebe, D.S., & Shustov, B.M., 

2002, ApJ, 569, 792 

Myers, PC, 2005, ApJ, 623, 280 
Myers, P.C. & Benson, P.J., 1983, ApJ, 266, 309 
Natarajan, P. & Lynden-Beh, D., 1997, MNRAS, 286,268 
Oberg, K. L, van Broekhuizen, F., Eraser, H. J., Bisschop, 

S. E., van Dishoeck, E. F., & Schlemmer, S. 2005, ApJL, 

621, L33 

Oberg, K. I., van Dishoeck, E. F., & Linnartz, H. 2009, 
AA, 496, 281 

Ogino, S., Tomisaka, K. & Nakamura, F., 1999, PASJ, 51, 
637 

Onsager, L., 1931, Phys Rev, 37, 405 
Ossenkopf, V., & Henning, T. 1994, AA, 291, 943 
Pagani, L., Bacmann, A., Cabrit, S., Vastel, C, 2007, AA, 
467, 179 

Pension, M.V., 1969, MNRAS, 144, 425 

Prasad, S. S., & Tarafdar, S. P. 1983, ApJ, 267, 603 

Rawlings, J.M.C., Hartquist, T.W., Menten, CM., 

Wilhams, D.A., 1992, MNRAS, 255, 471 
Rawlings, J.M.C. & Yates, J.A., 2001, MNRAS, 326, 1423 
Redman, M.P., Keto, E., Rawlings, J.M.C., 2006, MNRAS, 

370, 1 

Roberts, J. F., Rawlings, J. M. C., Viti, S., & Wilhams, 

D. A. 2007, MNRAS, 976 
Schnee, S., Caselli, P., Goodman, A., Arce, H. G., 

Ballesteros-Paredes, J., & Kuchibhotla, K. 2007, ApJ, 

671, 1839 

Shematovich, V.I., Wiebe, Shustov, B.M., & Li, Z.-Y., 

2003, ApJ, 588, 894 

Shen, C. J., Greenberg, J. M., Schutte, W. A., & van 

Dishoeck, E. F. 2004, AA, 415, 203 
Shu, F.H., 1977, ApJ, 214, 488 

Schutte, W. A., & Greenberg, J. M. 1991, AA, 244, 190 
Sohn, J. Lee, C.W., Park, Y-S., Lee, H.M., Myers, P.C, 

Lee, Y., 2007, ApJ , 664 928 
Spitzer, L., 1978, Physical Processes in the Interstellar 

Medium, New York Wiley-Interscience 
Tafalla, M., Myers, P., Caselh, P., Walmsley, C, & Comito, 

C. 2002, Ap J, 569, 815 
Takahashi, J., & Williams, D. A. 2000, MNRAS, 314, 273 
Tsamis, Y.G., Rawlings, J.M.C, Yates, J.A., & Viti, S., 

2008, arXiv:0803.0519 
van Dishoek, E. F., Blake, G. A., Draine, B. T., & Lunine, 

J. I. 1993, Protostars and Planets HI, 163 
Van Weeren, R.J., Brinch, C, & Hogerheijde, M.R., 2009, 

AA, 497, 773 

Ward- Thompson, D., Scott, P., Hills, R., & Andre, P. 1994, 



© 0000 RAS, MNRAS 000, 000-000 



Supercritical Starless Cores 9 

MNRAS, 268, 274 
Ward- Thompson, D., Motte, F., & Andre, P. 1999, MN- 
RAS, 305, 143 

Whitworth, A.R & Ward-Thompson, D., 2001, ApJ, 547, 
317 

Willacy, K., & Millar, T. J. 1998, MNRAS, 298, 562 
Williams, J. P., Myers, P.O., Wilner, D.J., DiEYancesco, J., 

1999, ApJ , 513, L61 
Williams, J.P., Lee, C.W., Myers, P.C., 2006, ApJ , 636, 

952 

This paper has been typeset from a TtiK/ M^jX file prepared 
by the author. 



© 0000 RAS, MNRAS 000, 000-000 



10 Keto and Caselli 




2.0 2.5 3.0 3.5 

log radius (AU) 



Figure 1. The progression of the density profile of a thermally 
supercritical starless core during subsonic contraction. The pro- 
file with a central density of 2 X 10^ cm~^ (the central density of 
L1544) is indicated by the dashed red line. The solid green line 
shows how well the approximation of equation [T] fits this density 
profile even though the approximation assumes isothermal gas 
and the numerical calculation allows for temperature variations 
due to radiative equilibrium (figure |3}. Density profiles at earlier 
and later times are shown in blue below and above this line. The 
evolution covers 1 Myr. For each profile, a cross marks the char- 
acteristic radius, rf (equation [2]l that separates the inner region 
with approximately constant density, and the outer region with 
density scaling as r~^. A slope of —2 is shown as the dashed green 
line. Over time, as the central density increases, the characteristic 
radius moves inward and the BE sphere resembles more and more 
a singular isothermal sphere with p ~ everywhere. 
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Figure 2. The progression of the velocity profile of a thermally 
supercritical core during subsonic contraction. The profiles have 
a characteristic A-shape with the maximum inward velocity near 
the characteristic radius, rf (equation |2]l . The profile at the time 
that the central density is 2 X 10^ cm~^ (the central density of 
L1544) is indicated by the dashed red line. 



o 
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Figure 3. The temperature, density, abundance, and inward ve- 
locity profiles of a contracting starless core at 3 times in the evo- 
lution of the contraction, when the central density is 2 X 10^ cm"'' 
(top), 2 X lO'^ cm~3 (middle), and 2 X 10* cm"^ (bottom). Red 
indicates the temperature in K, blue the log of the density in 
cm"'', green the ratio of the CO abundance to the undepleted 
abundance multiplied by 10. A ratio of 10 corresponds to no de- 
pletion and zero to total depletion. Purple indicates the inward 
velocity. In this figure a positive velocity is inward, and the veloc- 
ity units are 0.1 kms~^. The dashed lines show the temperature, 
density, and abundance of a static sphere with the same central 
density. The subsonically contracting sphere closely resembles the 
static sphere except for the velocities. 
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Figure 4. The abundance of CO for different rates of desorption. 
The dashed blaclc line shows the log of the density (cm~^). The 
other 4 curves show the abundance calculated from an equilib- 
rium between freeze-out, desorption, and photodissociation for 4 
different rates of desorption ( i|3.2l l. The lowest rate of desorption 
(blue curve) results in a CO abundance that is too low to match 
the observed CO line brightness. The highest rate of desorption 
(black curve) results in spectral line models that best match the 
observations. 



1.5 



E 



1.0 



0.5 



0.0 



-0.06 -0.04 -0.02 0.00 0.02 0.04 0.06 
Distance (PC) 

Figure 5. Observations of the trace of the integrated intensity of 
0^*0(1-0) (boxes and dotted line), C"O(l-0) (stars and dashed 
line), and N2H+(l- 0) (diamonds and solid line) across L1544 from 
ICaselli et al. 1 1I1999I ). The N2H+ molecule is not much depleted, if 
at all, at high densities, and its emission peaks toward the center 
of L1544. In contrast, the lines of CO, a molecule that is subject 
to depletion, are approximately constant across the core, possibly 
with lower intensity toward the center. The N2H+ emission has 
been divid ed by 5.65 and the C^*0 emission by 3.65. to match 
figure 2 in lCaselli et al.l l ll999l) 
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Figure 6. Trace of the integrated intensity of C^*O(l-0) (dotted 
line), C^^O(l-O) {dashed line), and N2H+(l-0) {solid line), across 
our model of a subsonically contracting core. This figure shows 
the integrated intensities at 3 times in the evolution of the con- 
traction, when the central density is 2 X 10^ cm""^ {top), 2 X lO'^ 
cm~'^ {middle), and 2 X 10* cm~^ {bottom). The internal struc- 
ture of this model is shown in figure |3] The N2H+ emission has 
been divided by 5.65 and the C^^^O emission by 3.65, the same as 
in figure [5] 
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Figure 7. Trace of the peak brightness of C^^O(l-O) {dotted line), 
0^^0(1-0) {dashed line), and N2H+(l-0) {solid line), across our 
model of a subsonically contracting core. This figure shows the 
peak brightness at 3 times in the evolution of the contraction, 
when the central density is 2 X lO'' cm~^ {top), 2 X 10*^ cm^"^ 
{middle), and 2 X 10* cm~'^ {bottom). The internal structure of 
this model is shown in figure |3] The N2H+ emission has been 
divided by 5.65 and the C^^O emission by 3.65, the same as in 
figure[5] The ratio of the peak intensity of 0^*0(1-0) and C'^''0{1- 
0) is different than the ratios of the integrated intensities because 
the C^^O(l-O) line is split into 3 hyperfine lines (figure |9j- 
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Figure 8. Trace of the FWHM of Ci*O(l-0) {dotted line), 
C^^O(l-O) {dashed line), and N2H+(l-0) {solid line), across our 
model of a subsonically contracting core. This figure shows the 
line width at 3 times in the evolution of the contraction, when 
the central density is 2 X 10^ cm"^ {top), 2 X 10^ cm^^ {middle), 
and 2 X 10* cm~'^ {bottom). The figure shows that in the cen- 
ter of the core the predicted C^*0, C^^O, and N2H+ line widths 
increase by 0.074, 0.060, and 0.145 kms~^, respectively because 
the unresolved infall velocities increase as the core evolves. The 
internal structure of this model is shown in figure \3\ 
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Figure 9. Spectra of CO(l-O) and N2H+(l-0) at the center of our model of a subsonically contracting core. This figure shows the line 
emission at 3 times in the evolution of the contraction, when the central density is 2 X 10^ cm~^, 2 X lO'^ cm~^, 2 X 10* cm~^. The upper 
left figure shows spectra of C^* O(l-0) (red), C^^O( l-O) (blue) at the 3 evolutionary times corresponding to the different densities. Also 
shown are the observed spectra l lCaseUi et al.lll999^ . The C^'^O spectra have been shifted upward by 1 K. There is very little difference 
in the model CO spectra as the core evolves. A slight decrease in peak intensity is noticeable from 2 X 10^ to 2 X 10* yr as the line 
width increases with the increasing infall velocities. The other 3 figures show the spectra of N2H"'"(l-0) at the 3 times with one figure for 
each time: 2 X 10^ cm~^ {upper right), 2 X 10^ cm"-^ (lower left), and 2 X 10* cm~^ (lower right). In these 3 figures, each N2H+ model 
spectrum (red) is shown together with the observed spectrum (blue). There is very little difference between the model spectra at the last 
two evolutionary times because the density increase is confined to a small region in the very center of the core. The internal structure of 
these models are shown in figure [3] 
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